function ax = glm_PlotBetas(hfile, c, plotopts)
% GLM::PlotBetas  - plot betas
%
% FORMAT:       [ax] = glm.PlotBetas(c [, plotopts])
%
% Input fields:
%
%       c           Cx3 coordinates to sample
%       plotopts    structure with optional fields
%        .ax        plot into given axis object
%        .cond      1xC condition selection (double or boolean)
%        .condcol   condition color (overrides predictor colors)
%        .cons      contrasts 1xC cell array with weights (from cond)
%        .csys      coordinate system (default: tal, see bvcoordconv)
%        .group     1xS double or 1xG cell array with group assignment
%        .legend    1xC cell array with (selected) condition names
%        .plottype  either of {'bar'}, 'scatter'
%        .regress   if given, either 'ols', or 'robust'
%        .scatter   values to scatter against (e.g. age, etc.)
%        .scrange   scatter range [x1, y1, x2, y2]
%        .shape     if only one coordinate given, shape sample around
%                   it, either of 'cube', 'sphere', {'voxel'}
%        .size      shape size, either cube edge length or radius
%        .se        plot standard error bars (if given and true)
%        .vmr       VMR object to link (for browsing)
%
% Output fields:
%
%       ax          axis object

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'glm') || ...
   ~isa(c, 'double') || ...
    isempty(c) || ...
    size(c, 2) ~= 3 || ...
    any(isinf(c(:)) | isnan(c(:)))
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end

% currently only VTC-based RFX GLMs supported
bc = bvqxfile_getcont(hfile.L);
if bc.ProjectType ~= 1 || ...
    bc.ProjectTypeRFX ~= 1
    error( ...
        'BVQXfile:BadArgument', ...
        'Only VTC-based RFX-GLMs supported for now.' ...
    );
end
numsubs = bc.NrOfSubjects;
numspred = bc.NrOfSubjectPredictors;

% get plot handles
ph = aft_Handles(hfile);
if ~isfield(ph, 'BetaPlot') || ...
   ~isstruct(ph.BetaPlot) || ...
   ~isfield(ph.BetaPlot, 'Axes') || ...
    numel(ph.BetaPlot.Axes) ~= 1 || ...
   ~ishandle(ph.BetaPlot.Axes)
    aft_SetHandle(hfile, 'BetaPlot', struct( ...
        'Axes',     [], ...
        'Bars',     [], ...
        'BarPos',   [], ...
        'Layout',   [0, 0], ...
        'Scatters', [], ...
        'ScLines',  [], ...
        'SEVLines', [], ...
        'SEULines', [], ...
        'SELLines', [], ...
        'Style',    'bar', ...
        'Text',     []));
    ph = aft_Handles(hfile);
end
ph = ph.BetaPlot;

% check options
if nargin < 3 || ...
   ~isstruct(plotopts) || ...
    isempty(plotopts)
    plotopts = struct;
else
    plotopts = plotopts(1);
end
if ~isfield(plotopts, 'ax') || ...
   ~isa(plotopts.ax, 'double') || ...
    numel(plotopts.ax) ~= 1 || ...
   ~ishandle(plotopts.ax) || ...
   ~strcmpi(get(plotopts.ax, 'type'), 'axes')
    if isfield(ph, 'Axes') && ...
        numel(ph.Axes) == 1 && ...
        ishandle(ph.Axes)
        ax = ph.Axes;
    else
        figure;
        ax = axes;
        ph.Axes = ax;
        ph.Bars = [];
        ph.Scatters = [];
    end
else
    ax = plotopts.ax;
    if isempty(ph.Axes) || ...
        ph.Axes ~= ax
        ph.Axes = ax;
        ph.Bars = [];
        ph.Scatters = [];
    end
end
if ~isfield(plotopts, 'cond') || ...
   (~isa(plotopts.cond, 'double') && ...
    ~islogical(plotopts.cond)) || ...
    isempty(plotopts.cond)
    plotopts.cond = 1:(numspred - 1);
end
if ~isfield(plotopts, 'condcol') || ...
   ~isa(plotopts.condcol, 'double') || ...
    size(plotopts.condcol, 1) ~= numel(plotopts.cond) || ...
    size(plotopts.condcol, 2) ~= 3 || ...
    any(isinf(plotopts.condcol(:)) | isnan(plotopts.condcol(:)) | ...
        plotopts.condcol(:) < 0 | plotopts.condcol(:) > 255)
    plotopts.condcol = [];
end
if ~isfield(plotopts, 'cons') || ...
   ~iscell(plotopts.cons)
    plotopts.cons = {};
end
for cc = numel(plotopts.cons):-1:1
    if ~isa(plotopts.cons{cc}, 'double') || ...
        isempty(plotopts.cons{cc}) || ...
        numel(plotopts.cons{cc}) ~= length(plotopts.cons{cc}) || ...
        numel(plotopts.cons{cc}) > numspred || ...
        any(isinf(plotopts.cons{cc}) | isnan(plotopts.cons{cc}))
        plotopts.cons(cc) = [];
        continue;
    end
    plotopts.cons{cc} = plotopts.cons{cc}(:)';
end
if ~isfield(plotopts, 'csys') || ...
   ~ischar(plotopts.csys) || ...
   ~any(strcmpi(plotopts.csys(:)', {'bvc', 'bvi', 'bvs', 'bvx', 'mni', 'tal'}))
    plotopts.csys = 'tal';
else
    plotopts.csys = lower(plotopts.csys(:)');
end
if strcmp(plotopts.csys, 'mni')
    c = round(mni2tal(c));
    plotopts.csys = 'tal';
end
if ~isfield(plotopts, 'group') || ...
   (~iscell(plotopts.group) && ...
    ~isa(plotopts.group, 'double'))
    plotopts.group = 'none';
end
scount = numsubs;
if isa(plotopts.group, 'double') && ...
   ~isempty(plotopts.group) && ...
    numel(plotopts.group) == max(size(plotopts.group)) && ...
   ~any(isinf(plotopts.group) | isnan(plotopts.group) | ...
        plotopts.group < 1 | plotopts.group > numsubs)
    plotopts.group = round(plotopts.group(:)');
    numgrps = max(plotopts.group);
    group = cell(1, numgrps);
    for gc = numgrps:-1:1
        group{gc} = find(plotopts.group == gc);
        if isempty(group{gc})
            group(gc) = [];
            continue;
        end
    end
    plotopts.group = group;
end
if iscell(plotopts.group) && ...
    numel(plotopts.group) < numsubs
    plotopts.group = plotopts.group(:)';
    scount = 0;
    for gc = numel(plotopts.group):-1:1
        plotopts.group{gc} = round(plotopts.group{gc}(:)');
        plotopts.group{gc}(isinf(plotopts.group{gc}) | isnan(plotopts.group{gc}) | ...
            plotopts.group{gc} < 1 | plotopts.group{gc} > numsubs) = [];
        if isempty(plotopts.group{gc})
            plotopts.group(gc) = [];
            continue;
        end
        scount = scount + numel(plotopts.group{gc});
    end
else
    plotopts.group = [];
end
if ~isfield(plotopts, 'legend') || ...
   ~iscell(plotopts.legend) || ...
    numel(plotopts.legend) ~= numel(plotopts.cond)
    plegend = cell(1, numel(plotopts.cond));
    for lc = 1:numel(plegend)
        plegend{lc} = regexprep(bc.Predictor(plotopts.cond(lc)).Name2, ...
            '^.*\:\s*(.*)\s*$', '$1');
    end
    plotopts.legend = plegend;
end
if ~isfield(plotopts, 'plottype') || ...
   ~ischar(plotopts.plottype) || ...
   ~any(strcmpi(plotopts.plottype(:)', {'bar', 'scatter'}))
    plotopts.plottype = 'bar';
else
    plotopts.plottype = lower(plotopts.plottype(:)');
end
if ~isfield(plotopts, 'regress') || ...
   ~ischar(plotopts.regress) || ...
   ~any(strcmpi(plotopts.regress(:)', {'ols', 'robust'}))
    plotopts.regress = '';
else
    plotopts.regress = lower(plotopts.regress(:)');
end
if ~isfield(plotopts, 'scatter') || ...
   ~isa(plotopts.scatter, 'double') || ...
    isempty(plotopts.scatter) || ...
    numel(plotopts.scatter) ~= length(plotopts.scatter) || ...
    numel(plotopts.scatter) ~= scount || ...
    any(isinf(plotopts.scatter) | isnan(plotopts.scatter))
    plotopts.scatter = [];
else
    plotopts.scatter = plotopts.scatter(:);
end
if ~isfield(plotopts, 'scrange') || ...
   ~isa(plotopts.scrange, 'double') || ...
    numel(plotopts.scrange) ~= 4
    plotopts.scrange = [-Inf, Inf, -Inf, Inf];
else
    plotopts.scrange = plotopts.scrange(:)';
end
if isfield(plotopts, 'se')
    drawse = true;
else
    drawse = false;
end
if isfield(plotopts, 'vmr') && ...
    numel(plotopts.vmr) == 1 && ...
    isBVQXfile(plotopts.vmr, 'vmr')
    linkvmr = true;
else
    linkvmr = false;
end

%msz = size(bc.GLMData.RFXGlobalMap);
mnl = numel(bc.GLMData.RFXGlobalMap);

% get coordinate indices
if ~strcmp(plotopts.csys, 'bvx')
    c = bvcoordconv(c, [plotopts.csys '2bvx'], aft_BoundingBox(hfile));
end

% get map contents
bv = zeros(numsubs, numel(plotopts.cond));
for sc = 1:numsubs
    for mc = 1:numel(plotopts.cond)
        bv(sc, mc) = mean(bc.GLMData.Subject(sc).BetaMaps( ...
            c + (plotopts.cond(mc) - 1) * mnl));
    end
end

% build contrasts
if ~isempty(plotopts.cons)
    cbv = zeros(numsubs, numel(plotopts.cons));
    for cc = 1:numel(plotopts.cons)
        nci = min(numel(plotopts.cons{cc}), size(bv, 2));
        cbv(:, cc) = sum(repmat(plotopts.cons{cc}(1, 1:nci), [numsubs, 1]) .* ...
            bv(:, 1:nci), 2);
    end
    bv = cbv;
end

% group subject data?
if iscell(plotopts.group)
    bvm = zeros(numel(plotopts.group), size(bv, 2));
    bvs = bvm;
    for gc = 1:numel(plotopts.group)
        bvm(gc, :) = mean(bv(plotopts.group{gc}, :), 1);
        if drawse
            try
                bvs(gc, :) = std(bv(plotopts.group{gc}, :), [], 1) ./ ...
                    sqrt(numel(plotopts.group{gc}));
            catch
                % do nothing
            end
            %disp(bv(plotopts.group{gc}, :));
        end
    end
else
    bvm = mean(bv, 1);
    if drawse
        try
            bvs = std(bv, [], 1) ./ sqrt(size(bv, 1));
        catch
            bvs = zeros(size(bvm));
        end
    end
end

% plot type
switch plotopts.plottype
    
    % bar plots
    case {'bar'}

        % plot either update bars or draw new
        nbg = size(bvm, 2);
        onlyone = false;
        if size(bvm, 1) == 1
            onlyone = true;
            bvm(2, :) = 0;
            bvs(2, :) = 0;
        end
        if isempty(ph.Bars)
            bh = bar(ax, bvm);
            ph.Bars = bh;
        else
            for rc = 1:size(bvm, 2)
                set(ph.Bars(rc), 'YData', bvm(:, rc)');
            end
        end

        % draw standard error bars
        if drawse
            bmu = bvm + bvs;
            bml = bvm - bvs;
            bps = [0.5, 0.6, 2/3, 0.72, 0.76, 0.78, 0.79, 0.8];
            bpi = bps(min(8, size(bvm, 2)));
            if isempty(ph.SEVLines)
                hold(ax, 'on');
                barpos = (1:size(bvm, 1))' * ones(1, nbg) + ...
                    ones(size(bvm, 1), 1) * ((-bpi / 2 + bpi / (2 * nbg)):bpi / nbg:bpi / 2);
                ph.BarPos = barpos;
                barpol = barpos - 0.2 * bpi / nbg;
                barpor = barpos + 0.2 * bpi / nbg;
                hold(ax, 'off');
                vln = zeros(size(barpos));
                uln = vln;
                lln = vln;
                for lc = 1:numel(barpos)
                    vln(lc) = line([barpos(lc), barpos(lc)], [bml(lc), bmu(lc)]);
                    uln(lc) = line([barpol(lc), barpor(lc)], [bmu(lc), bmu(lc)]);
                    lln(lc) = line([barpol(lc), barpor(lc)], [bml(lc), bml(lc)]);
                end
                ph.SEVLines = vln;
                ph.SEULines = uln;
                ph.SELLines = lln;
            else
                nbg = size(bvm, 2);
                barpos = ph.BarPos;
                barpol = barpos - 0.2 * 0.8 / nbg;
                barpor = barpos + 0.2 * 0.8 / nbg;
                vln = ph.SEVLines;
                uln = ph.SEULines;
                lln = ph.SELLines;
                for lc = 1:numel(barpos)
                    set(vln(lc), 'XData', [barpos(lc), barpos(lc)], 'YData', [bml(lc), bmu(lc)]);
                    set(uln(lc), 'XData', [barpol(lc), barpor(lc)], 'YData', [bmu(lc), bmu(lc)]);
                    set(lln(lc), 'XData', [barpol(lc), barpor(lc)], 'YData', [bml(lc), bml(lc)]);
                end
            end
            if onlyone
                set(ax, 'XLim', [0.5 1.5]);
            end
        end
        
        % scatter
    case {'scatter'}
        
        % plot either update bars or draw new
        if isempty(ph.Scatters)
            sh = scatter(ax, plotopts.scatter, bv);
            ph.Scatters = sh;
        else
            set(ph.Scatters, 'YData', bv);
        end
        
        % regression
        scm = mean(plotopts.scatter);
        scn = size(bv, 1);
        sX = [plotopts.scatter - scm, ones(scn, 1)];
        scrange = plotopts.scrange;
        if isinf(scrange(1))
            scrange(1) = min(plotopts.scatter);
        end
        if isinf(scrange(2))
            scrange(2) = min(bv);
        end
        if isinf(scrange(3))
            scrange(3) = max(plotopts.scatter);
        end
        if isinf(scrange(4))
            scrange(4) = max(bv);
        end
        switch plotopts.regress
            case {'ols'}
                regb = pinv(sX' * sX) * sX' * bv;
            case {'robust'}
                [regb, regr, sw, regs] = fitrobustbisquare(sX, bv);
                wc = get(ph.Scatters, 'Children');
                for wcc = 1:numel(wc)
                    wf = 1 - sw(end+1-wcc);
                    set(wc(wcc), 'MarkerEdgeColor', [wf, 0, 1 - wf]);
                end
                R2 = max(0, 1 - var(sw .* regr) / var(sw .* bv));
                R2s = sprintf('R2 = %.3f\np = %.4f', R2, regs.p(1));
                if regb(1) > 0
                    py = 0.15;
                else
                    py = 0.85;
                end
                px = scrange(1) + 0.85 * (scrange(3) - scrange(1));
                py = scrange(2) + py * (scrange(4) - scrange(2));
                if isempty(ph.Text)
                    tx = text(px, py, R2s);
                    ph.Text = tx;
                else
                    set(ph.Text, 'String', R2s, 'Position', [px, py, 0]);
                end
            otherwise
                regb = [];
        end
        if ~isempty(plotopts.regress)
            ly1 = regb(2) + (scrange(1) - scm) * regb(1);
            ly2 = regb(2) + (scrange(3) - scm) * regb(1);
            if isempty(ph.ScLines)
                lh = line(scrange([1, 3])', [ly1; ly2]);
                ph.ScLines = lh;
            else
                set(ph.ScLines, 'XData', scrange([1, 3]));
                set(ph.ScLines, 'YData', [ly1, ly2]);
            end
        end
        set(ax, 'XLim', scrange([1, 3]), 'YLim', scrange([2, 4]));
end

% update UI
drawnow;

% set handles
aft_SetHandle(hfile, 'BetaPlot', ph);

% link to VMR
if linkvmr
    vmr = plotopts.vmr;
    plotopts = rmfield(plotopts, 'vmr');
    plotopts.ax = ax;
    vh = aft_Handles(vmr);
    vh.PlotHandles.BetaPlotGLM = hfile;
    vh.PlotHandles.BetaPlotOpts = plotopts;
    aft_SetHandle(vmr, 'PlotHandles', vh.PlotHandles);
    aft_Browse(vmr);
end
